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r~| Abstract 

i A recently developed method for the calculation of Lyapunov exponents of dynam- 

q ical systems is described. The method is applicable whenever the linearized dynamics 

is Hamiltonian. By utilizing the exponential representation of symplectic matrices, 
this approach avoids the renormalization and reorthogonalization procedures neces- 
sary in usual techniques. It is also easily extendible to damped systems. The method 
is illustrated by considering two examples of physical interest: a model system that 
describes the beam halo in charged particle beams and the driven van der Pol oscil- 
lator. 



e-mail: 

habib@predator.lanl.gov 

ryne@lanl.gov 



1 Introduction 

Chaotic dynamical systems exhibit exponential divergence of initially nearby trajecto- 
ries. This divergence is quantified by the Lyapunov exponents of the system which are 
obtained from linearizing the dynamics around a fiducial trajectory [1]. Over the past 
two decades or so there has been "intense activity" [2] directed toward the compu- 
tation of these exponents resulting in several different numerical approaches [1] [3] [4] . 
The two obvious difficulties associated with the computation of Lyapunov exponents 
are: (1) exponential growth of the separation vector (between the fiducial and nearby 
trajectories) and (2) the exponential collapse of initially orthogonal separation vectors 
onto the direction of maximal growth. Most conventional methods overcome these 
hurdles by intermittent numerical rescaling and reorthogonalization (through, e.g., 
the Gramm-Schmidt procedure [3]). Many chaotic systems are Hamiltonian or they 
can be transformed into a Hamiltonian system by suitable manipulations. However, 
none of the above general methods are designed to take advantage of this fact. 

As is well known, the dynamics of classical Hamiltonian systems has an underly- 
ing symplectic structure [5]. In recent years symplectic methods have been applied 
with great success to classical dynamical problems. The field of accelerator dynamics 
has been revolutionized by the introduction of nonlinear symplectic maps as repre- 
sented by Lie transformations [6] [7]. Very long time integration of charged particle 
and planetary systems has been aided by the development of high order symplectic 
integration algorithms [8]. Recently a symplectic map-based approach for the calcu- 
lation of Lyapunov exponents has been developed [9]. As shown below, this approach 
obviates analytically the need for rescaling and reorthogonalization in the numerical 
computation of the exponents. 



2 The Method 

Consider a 2-m dimensional continuous-time dynamical system governed by the equa- 
tions 

£ = F(.,0, (1) 

where z = (zi, z 2 , ■ ■ ■ , £ 2m ) an d similarly for F. Let z denote some given fiducial 
trajectory. Define deviations from this trajectory by letting Z = z — z , and linearize 
the above equations. The new set of equations for the deviation variables is 

^ = DF(z ,t)-Z. (2) 
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The approach described below can be used whenever this linearized set of equations 
is derivable from a Hamiltonian. From now on suppose that this is the case, and that 
one can write Z = (gi, q 2 ,- ■ ■ , q m ,Pi,P2, ' ' ' iVm), where qi and pi denote canonically 
conjugate coordinates and momenta, respectively. It follows that 

f = -WZ}, (3) 

where {,} denotes the Poisson bracket, and where if is a homogeneous quadratic 
polynomial in the qi and pi. A system such as this is governed by a symplectic matrix 
M that maps the initial variables Z m into time-evolved variables Z(t), 

Z(t) = M(t)Z m (4) 

Let A be given by 

A = lim (MM) , (5) 

where M denotes the matrix transpose of M. The Lyapunov exponents then equal 
the logarithm of the eigenvalues of A [1]. 

It is easy to show that M satisfies the equation of motion (See, for example, Ref. 

[11]) 

where S denotes the symmetric matrix given by 

H{Z,t) = -Y l S ii Z i Z i , (7) 

and where 

Here 1 denotes the m X m identity matrix. It follows that the evolution of MM is 
governed by the equation 

—MM = JSMM - MMSJ. (9) 

dt y ' 

Standard methods for obtaining the Lyapunov exponents deal with M, which is not 
real symmetric (hence the need for reorthogonalization) and which has exponentially 
growing elements. To avoid these difficulties one exploits the fact that M is symplectic 
by making use of the exponential representation of symplectic matrices [6]: Any 
symplectic matrix M can be written in the form 

M = e JSa e JS < (10) 



where S a is a symmetric matrix that anticommutes with J and S c is another sym- 
metric matrix that commutes with J. It is important to note that the second matrix 
on the right hand side of (10) is in fact unitary, so that 

MM = e 2JS \ (11) 

Note that this matrix has fewer degrees of freedom than M and its eigenvectors are 
orthogonal. Rather than attempting to directly integrate Eqn. (9) which would still 
have a large numbers problem, focus attention instead on the exponent JS a in Eqn. 
(10). It is clear that now there is no large numbers problem since JS a already appears 
as an exponent. 

To proceed further, one obvious approach is to use an explicit representation of 
exp( JS a ). Such a representation is well known for Sp(2) and has recently been found 
for Sp(4:) (see the Appendix). Generalizations to Sp(2m) are in progress [10]. To 
illustrate the method it is convenient to restrict attention to systems with a two- 
dimensional phase space. When driven, these represent the simplest continuous-time 
systems that can exhibit chaos. The most general two-dimensional symplectic matrix 
can be written in the form 

M = e JSa e JSc 

__ /J.(B2Cosa+B 3 s'm.a) bBi (1 2^ 



where a, b and \i are real coefficients and where the B{ are basis elements of the Lie 
algebra sp{2) [6]: 

Bi = I \ : ) , b 2 = ( ° I ) , b 3 = I : \ ] . (i3) 
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It follows that 



Thus, one obtains, 



MM = e 2 M- B 2Cosa+B 3 sina) Q^N 



A = lim e 0V*)( B 2 cos a-HB 3 sin a) Qg\ 

t-K» ' ^ ' 

Finally, it is easily shown that the eigenvalues of this matrix are e ±M ' t . The Lyapunov 
exponents are then equal to ±///t in the limit t — > oo. With this convenient choice 
of variables, the explicit representation of MM is given by 

(cosh lu, + sin a sinh 2u cos a sinh 2u \ , 

. . • (16) 

cos a sinh 2/x cosh 2/x — sin a sinh 2/x / 

The unknown quantities a and \i can grow in time at most as 0(t). Differential 
equations for these quantities can be obtained by returning to Eqn. (9), the dynamical 
equation for MM. 
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For simplicity, assume that H contains no term proportional to qp, so that the 
matrix S in (7) is of the form 

(17) 




After some manipulation, Eqns. (7)-(9) lead to the following: 

d\i 1 

~dl 
da 



, -^22 - 511 J cos a, 



: 5n + 5 2 2 — (^22 — ■Sn)sinacoth/x. (18) 

From the initial condition M(0) = /, if one chooses /i(0) = 0, then cos 2 a(0) = 1, i.e., 
a(0) = or 7r. These differential equations form the basis of the method for calculating 
the Lyapunov exponents of Hamiltonian systems: They are stepped forward in time 
numerically till some desired convergence for the exponents, ±///t, is achieved. It will 
be shown later how to apply the method to certain non-Hamiltonian systems. 

3 Applications: Two Examples 

As a first concrete example, consider the newly developed "core-halo" model which 
describes beam halo in mismatched charged particle beams [12]. The transverse 
equation of motion for a halo particle in this model, assuming constant external 
focusing, is 

x + x-{l- V 2 )f{x,r{t)) = (19) 

where x is the position variable for a halo particle, f(x,r(t)) is the force due to the 
space charge of the beam core, and r(t) is the time dependent rms radius of the core. 
The core radius is assumed to follow the envelope equation 

r + r-l--?--'£ = 0. (20) 

Here units have been chosen so that the time independent solution of (20) (i.e., a 
matched beam) is given by r = 1. In these units rj = corresponds to the space charge 
dominated regime, while rj = 1 corresponds to the emittance dominated regime. Now 
assume 

/M = ^ (21) 

which has the correct asymptotic behavior: the force is linear when x <C r, and it is 
inversely proportional to x when x ^> r. The Eqns. (19) and (20) describe a driven 
nonlinear system with a mixed phase space as demonstrated by the stroboscopic 
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Figure 1: Stroboscopic plot of the chaotic sea in the core-halo model. Snapshots were 
taken at successive beam minima for 32 test particles. Parameter values were r(0) = 0.6, 
r(0) = 0, and 77 = 0.2. 
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Figure 2: Positive Lyapunov exponent for the core-halo model in a typical run. Parameters 
are the same as in Fig. 1. The simulation was run for 10 5 periods of the driving force 
using 100 integration steps per period for Eqns. (18) with a third-order Runge-Kutta 
algorithm. The envelope equation and the fiducial trajectory were integrated with a fourth- 
order symplectic algorithm using 200 steps per period. 



plot shown in Fig. 1. The presence of a chaotic band is important because particles 
initially in the core can leak through the broken separatrix and be carried to large am- 
plitudes. The presence of such large amplitude particles can cause unacceptably high 
radioactivation levels in high intensity linacs planned for future accelerator-driven 
technologies [13]. Leakage through the separatrix can be enhanced through particle 
collisions and recent work has shown that this rate is controlled by the Lyapunov ex- 
ponent [14]. The Lyapunov exponent for this system may be computed by integrating 
(18) with 

3ii = l-(l-i/ 2 ^ l 2XI 



s 22 = 1 (22) 

where x denotes the fiducial trajectory. Fig. 2 displays the result for the Lyapunov 
exponent against time. The slow convergence of the exponent to its asymptotic value 
is typical of Hamiltonian systems. 

So far only explicitly Hamiltonian systems were considered. However, the only 
real requirement for using the method is that the linearized deviation equations in 
some variables be Hamiltonian. This allows for the inclusion of damped systems in 
the scheme. As an example, consider the following general driven nonlinear oscillator 



x + \(l- ex 2 ) x + V'(x) = a cos(tut). (23) 

By appropriate choices of A, e, and V(x), this reduces to an assortment of well- 
known equations including van der Pol (A < 0, e = 1, V(x) = (l/2)x 2 ), Duffing 
(A > 0, e = 0, V{x) = ax 2 + /3x 4 ), and the damped driven pendulum (A > 0, e = 0, 
V(x) = 1 — cos(x)). In terms of the deviation variable 8, the linearization of (23) 
yields 

6 + A (l - ex 2 ) 6 + (V"(x ) - 2e\x x ) 6 = (24) 



where x represents the fiducial trajectory. Introducing the new variable A defined 
through 

A = 6e~ 9 ^ (25) 



where 



4 A - e <) , (26) 



Eqn. (24) reduces to that describing an undamped oscillator with time dependent 
frequency, 

A + (V"(x ) - e\x x - 1 A 2 (l - ex 2 ) 2 ) A = 0. (27) 



It is now straightforward to proceed in the usual way: for linear damping (e = 0) the 
Lyapunov exponents %± of this system are given by 

X± = -^A±lim^ (t) (28) 

l t->oo t 

where /x follows from solving (18) for the system defined by (27). When e ^ 0, 

X ± = lim - (g(t) ± /x (0) (29) 

t->oo t 

modulo terms that are exponentially suppressed at late times. 

Figs. 3(a) and 3(b) show the Lyapunov exponents of the van der Pol system. For 
the chosen set of parameters these results are in agreement with those of Ref. [4]. In 
contrast with the results shown in Fig. 2, the convergence of the exponents is much 
faster, as is typical of non-weakly damped systems. 



4 Conclusion 

To summarize, a method for computing Lyapunov exponents that exploits the un- 
derlying symplectic structure of Hamiltonian dynamics has been reviewed. Just as 
symplectic integrators are not a panacea for all time integration problems, this method 
does not have universal applicability and advantages. However, when applicable, the 
method has certain advantages over standard techniques, most importantly the lack 
of systematic errors associated with intermittent reorthogonalization and rescaling. 
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Figure 3: (a)Positive Lyapunov exponent for the van der Pol oscillator with parameters 
A = —5, a = 5, and u = 2.466 (taken from Ref. [4]). The simulation was run for 10 5 periods 
of the driving force using 100 million time steps for Eqns. (18) and 200 million time steps 
for the fiducial trajectory. The integrators were third-order Runge-Kutta and fourth-order 
symplectic, respectively. 
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Figure 3: (b)Negative Lyapunov exponent for the van der Pol oscillator with the same 
parameters as in Fig. 3(a). 
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A Appendix: Summary of Sp(4) Results 

There exists a particularly convenient way of organizing the basis elements of sp(4) 
in terms of a triplet of matrices (the F{, Gi, and the B{) each triplet consisting in 
turn of three matrices. The B{ belong to the unitary sector which is irrelevant to the 
computation of the Lyapunov exponents. The most general (4 X 4) symplectic matrix 
then turns out to be of the form exp(a • F + b • G) times a unitary matrix, where 

a = {ai,a2,a 3 }, b = {&i,& 2) & 3 } 

F = {F 1 ,F 2 ,F 3 }, G = {G 1 ,G 2 ,G 3 }. (30) 

Here, the a^ and the b{ are scalars, and the explicit forms of the matrices F{ and G{ 
are given below. 

It turns out to be possible to resum the formal exponential above and to obtain 
the following result [10]: 

M = exp(a-F + b-G) 

= -(coshf^] + cosh.[u y ])I 

+ -(sch[u x ] + sch[u y ]) (a • F + b • G) 

4 
+ r~2 jr( cos h[w x ] - cosh[it v ]) (a x b) • K 

\ u x ~ U y) 

+ T~? ^(schH - sch[u y ]) ([b x (a x b)] • F + [a x (b x a)] • G) (31) 

\ U l ~ U y) 

where "•" represents the usual vector dot product, "x" represents the vector cross 
product, I is the (4 X 4) identity matrix and 

sch[x] = sinh[x]/x, 
u x = V2s 2 + 4v, 

Uy = \J1S 2 — 4f, 



s = a • a + b • b, 

v 2 = (axb)-(axb). (32) 

A new triplet of matrices, the K{ appears. This is because the F{ and the G{ do not 
form a closed subalgebra. The K{ are idempotent and unitary. Note that all the F{, 
Gi and Ki are traceless. Thus, quite trivially, Tr(M) = 2(cosh[/u x ] + cosh[/u y ]) and 
the eigenvalues of M are exp[/u x ], exp[— u x ], exp[it v ], exp[— u y ]. 
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Finally, the explicit forms of the above matrices are given by: 



Ki 



( -1 \ 

10 

0-10 

1 



,K 2 



( -1 \ 
-10 
0-1 
0-10 



,K 3 



( -1 \ 

10 

10 

-10 



(33) 



Fi 



/ -1 -i \ 

-10-10 

0-101 

-10 10 



,F 2 



/ 1 1 \ 
0-10-1 
10-10 
0-101 



,*3 



/ 1 -1 \ 

10-1 

-10-10 

0-10-1 

(34) 



/ 



Gi 







-1 1 \ 
-10 10 

10 1 

1 1 j 



,G 2 



( 1 -1 \ 

0-101 

-10-10 

10 1 



,G 3 



( -1 -1 \ 
0-10-1 
-10 10 
0-101 



(35) 
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